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ABSTRACT 


A  numerical  solution  of  the  laminar,  steady  near  wake 
of  an  axisymmetric  body  in  supersonic  flow  is  developed. 

The  body  considered  is  a  spherically  capped  cylinder  with 
a  truncated  base.  The  full  set  of  compressible  partial 
differential  equations  are  solved. 

The  numerical  approach  is  a  two-step  time  dependent 
finite  difference  procedure  similar  to  that  used  by  Chang 
and  Allen  in  two-dimensions.  The  primary  feature  of  this 
scheme  is  the  elimination  of  the  viscosity  dependent  sta¬ 
bility  condition  for  the  time  dependent  solution.  A 
variable  mesh  grid  system  is  included  to  provide  higher 
resolution  in  areas  of  flow  nonuniformity  and  both  the 
equation  system  and  the  finite  differencing  are  done  in 
conservation  form.  The  difference  forms  at  the  boundaries 
are  obtained  by  means  of  telescoping  with  the  internal 
difference  forms. 

A  series  of  numerical  examples  were  calculated  on  a 
CDC  6600  computer  with  both  a  coarse  and  fine  mesh  config¬ 
uration  for  various  free  stream  Reynolds;  numbers. 

Convergent,  steady  state  solutions  were  obtained  which 
exhibited  definite  recirculation  regions  and  showed  evidence 
of  recompresaion  shocks.  Restrictions  on  the  Reynolds 
number  for  a  given  mesh  size  were  investigated,  as  was  the 
relationship  between  solution  curvatures  and  mesh  size. 
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NOMENCLATURE 


a  General  coefficients  in  difference  forms  and  constant 
in  Sutherlands  Law 

d  Constant  in  isentropic  pressure  density  relationship 
e  Energy  normalized  by  free  stream  value 

e  Normalized  energy  resulting  from  first  iteration  step 

f  General  function  in  difference  forms 


Fx  , 


Fa, 
F.  , 
F* 


► 

i 


Functions  appearing  on  right  hand  side  of  difference 
form  of  continuity,  axial  momentum, radial  momentum 
and  energy  equations  (Equations  8  and  9). 


F  Portions  of  F  functions  consisting  of  non-central  parts 
of  twice  spatially  differenced  diffusive  terms 


Fa , 
Fa, 
F* 


F  functions  for  continuity,  axial  momentum,  radial 
momentum  and  energy  equations 


G  Functions  containing  updated  central  portions  of  twice 
spatially  differenced  diffusive  terms 

Ga,l  G  functions  appearing  in  denominators  of  difference 

.  I  form  of  axial  momentum,  radial  momentum  and 
s'f  energy  equations. 

Gi  J 

H  Portions  of  F  functions  which  are  not  twice  spatially 
differenced  diffusive  terms  and  constant  stagnation 
enthalpy  along  streamline 

k  Thermal  conductivity  normalized  by  free  stream  value 

k  Normalized  thermal  conductivity  resulting  from  first 

iteration  step 

M  Mach  number 
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90 


M  =  ua /e 

*  <*  • 

M  =  1/y  Ma 

m  'so 

P  Normalized  static  pressure 

Pr  Prandtl  number 

q  =  (uz+v2)Ma/2 

r  Radial  coordinate  normalized  by  base  half  height 

P.  Gas  constant 

Re  Reynolds  number  based  on  base  half  height 

Re^  Reynolds  number  based  on  mesh  size 
t  Time  normalized  by  base  half  height  uw 

u  Axial  velocity  component  normalized  by  free  stream  velocity 

u  Normalized  axial  velocity  resulting  from  first  iteration 
step 

v  Radial  velocity  component  normalized  by  free  stream 
velocity 

v  Normalized  radial  velocity  resulting  from  first  iteration 

step 

x  Axial  coordinate  normalized  by  base  half  height 
X  -  (y  -  l)/(ym  -  1; 

a  Ratio  of  adjacent  axial  mesh  sizes 

P  Ratio  of  adjacent  radial  mesh  sizes 

y  Ratio  of  specific  heats 

0  natio  of  mesh  sizes  Ar  and  Ar 

m+2  m 

Ar  Radial  mesh  size 
At  Normalized  time  step 
Ak  Axial  mesh  size 

e  Limiting  value  for  stability  criteria 
C  .  (y-l)R/(y  -1) R 

C  Value  of  C  resulting  from  first  iteration  step 
0  Flow  deflection  along  upper  boundary 
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4  Fluid  viscosity  normalized  by  free  stream  value 
a  Normalized  viscosity  resulting  from  first  iteration  step 
v  Kinematic  viscosity 

p  Density  normalized  by  free  stream  value 

p  Normalized  density  resulting  from  first  iteration  step 

Subscripts 

A  Starting  point  of  upper  boundary  characteristics 
extrapolation 

B  Intersection  point  of  first  family  characteristic  line 

and  upper  boundary  line 

c  Centerline  value 

C  Intersection  point  of  streamline  and  x  ■  xA  line 
m  Radial  mesh  point  index 

M  Upper  boundary  line  index 

n  Axial  mesh  point  index 

N  Downstream  boundary  line  index 

w  Body  wall  va lue 

•  Free  stream  value,  ahead  of  bow  shock 

0  Point  one  half  mesh  from  wall 

1  Point  one  and  one  half  meshes  from  wall  or 
point  one  mesh  from  f, 

2  Point  two  meshes  from  f, 

Superscript 
j  Time  step  index 
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INTRODUCTION 


This  work  aims  at  treating  the  laminar,*  steady,  axisymmetric 
near  wake  of  a  body  with  a  cylindrical  base  traveling  at  super¬ 
sonic  velocity  by  the  numerical  solution  of  the  full  set  of  com¬ 
pressible  partial  differential  equations  governing  the  motion. 

This  approach  eliminates  the  need  for  approximations,  such  as  the 
boundary  layer  approximation  and  the  patching  of  separate  regions 
where  different  approximations  are  made.  Numerical  solutions  also 
have  the  advantage  that  they  give  all  the  details  of  the  flow, 
while  methods  such  as  approximate  integral  solutions  do  not.  In 
theory  at  least,  numerical  solutions  can  be  used  to  check  their 
own  accuracy  by  making  calculations  at  successively  smaller  grid 
sizes. 

The  differential  equations  under  consideration  are  the  full 

compressible  Navier  Stokes  equations  in  axisymmetric  coordinates  and 

the  numerical  scheme  consists  of  a  two-step  time -dependent  finite 

difference  procedure.  Several  such  techniques  have  been  considered,1 

3 

but  the  most  promising  appears  to  be  the  approach  of  Cheng  and 
4  5  4 

Allen.  '  Allen  has  treated  the  two-dimensional  case,  utilizing  a 
uniform,  but  non-square  mesh  and  assuming  constant  viscosity  and  con¬ 
ductivity  and  simplified  outside  flow.  The  present  work  attempts  to 
extend  this  to  an  axisymmetric  configuration  and  realistic  outer  flow 
Moreover,  a  nonuniform  mesh  is  being  introduced  to  give  greater  reso¬ 
lution  in  the  high-gradient  regions. 

1,2 

The  other  possibly  competitive  techniques  have  been  consideret 
and  discarded.  Schemes  with  viscosity  dependent  stability  criteria, 
such  as  Reference  1,  appear  to  be  susceptible  to  numerical  diffi¬ 
culties  in  the  body  wall  and  base  regions,  those  difficulties  taking 
the  form  of  very  small  or  even  negative  density  values  as  the 
time -evolution  of  the  flow  computation  proceeds.  Such  phenomena 

*  The  feasibility  of  extending  this  work  to  turbulent  flows  is 
discussed  in  Appendix  III. 
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play  havoc  with  the  stability  of  the  computation.  The  scheme  of 

Reference  2  was  found  to  have  stability  difficulties  at  the  outflow 

4 

and  possibly  the  upper  boundaries.  We  are  attempting  to  take 
advantage  of  the  numerical  discovery  by  Cheng  and  Allen,  of  these 
pitfalls  and  of  the  techniques  which  appear  to  avoid  their 
occurrence . 


TR-763 


2 


MODEL  PROBLEM  AND 


INTERNAL  NUMERICAL  SCHEME 
Configuration 

The  type  of  axisymmetric  body  being  considered  is  shown  in 
Figure  1.  Note  that  the  base  section  is  cylindrical,  which  is 
a  simplifying  assumption  of  this  program.  The  calculation  domain 
being  considered  is  ABCFGHJA.  The  inflow  boundary  AC  is  set 
back  upstream  of  the  base  HJ  to  permit  upstream  influence  in  the 
boundary  layer  to  occur.  The  bow  shock  and  the  outer  region  CDEF 
are  computed  by  the  method  of  characteristics,  where  character¬ 
istic  BF  defines  the  limit  of  influence  of  the  boundary  layer 
up  to  point  A.  This  characteristics  calculation  is  used  to 
obtain  the  inviscid  portion  of  the  flow  field  upstream  of  the 

line  ABCD.  The  viscous  region  upstream  of  line  AB  is  obtained 

6 

from  a  boundary  layer  calculation,  which  is  started  far  upstream 
of  point  A  with  an  initial  step  profile.  The  downstream  boundary 
GF  is  located  where  the  flow  has  become  entirely  supersonic. 

The  data  on  this  boundary  is  obtained  by  an  extrapolation  from 
the  calculated  flow  field. 

Spatial  Differencing 

Since  a  steady  state  solution  is  the  objective,  it  is 
necessary  to  develop  a  consistent  set  of  spatial  differencing 
approximations  for  the  steady  state  terms.  The  flow  field  being 

considered  is  expected  to  contain  discontinuities  such  as  shocks. 

.  7 

Therefore  it  is  desirable  as  pointed  out  by  Lax,  Wendroff  et  al. 

to  use  the  conservation  form  of  the  differential  and  difference 
equations.  In  principle  the  conservation  form  of  the  Navier 
Stokes  equations  in  three-dimensional  cartesian  coordinates  can 
be  retained  under  a  transformation  to  three-dimensional  cylindric 
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coordinates.®  However,  when  this  conservation -preserving  set  of 
equations  is  reduced  by  the  assumption  of  axisymmetry,  circum¬ 
ferential  terms  and  derivatives  vanish  except  for  a  circumferential 
normal  stress  term  in  the  radial  momentum  equation.  This  term 
cannot  be  put  into  conservation  form  but  this  should  not  be 
disastrous  since  there  will  be  no  jumps  across  discontinuities 
in  the  circumferential  direction.  The  resulting  non-dimensional- 
ised,  unsteady  continuity,  axial  and  radial  momentum  and  energy 
equations  are: 


(pr)  t  +  (oru)x  +  (prv)r  * 


(pru) t  +  Qjruv  -  jjj-  (vx  +  +  £pr(ua  +  ft^eX) 


4  ru 

-  —  „  u 
3  Re  x 


2  ru  ,  ,  v., 

+  T  „  "  (v  +  — ) 

3  Re  r  r  jx 


=  0 


(prv)t  + 


pruv 


_  at_ 


_  (v  +  u  )  : 
Re  x  r  jx 


£r<v»+A.eX)  -f“-vr+f  ^-(ux^)]r 

OD  tO 


-  Pft-ex  -  7  *h(vr+ux)  +  3  it  r  ■  0 


(1) 
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Y 


[pr(e+q)]t  +  ^Dru(Ye+q)  -  frk(Ce)x] 

r  oo 

M 

-  rT  C  3  rtj(u,)x  +  ft(V,)x  ■  3  ruUVr+r,jVUr-  3  UUV]1 


Y 

+  -f  orv  (Ye+q)  -  p"—  [rk(^e)r3 

r  co  “ 


-iT[f  r“(v,)r+f“<ua>r 
00 


2 

J  ruvu 


+  ruuv 


-  3  yV 


where 


Y  -1 


q  *(u8  +  v3)  — 
M  -  u*  /e 

CD  mm 


(1)  concluded 


y  is  the  ratio  of  specific  heats,  R  is  the  gas  constant.  Re  the 
Reynolds  number,  Pr  the  Prandtl  number,  and  all  flow  variables  are 
normalized  by  thei?  values  in  the  free  stream,  which  are  denoted 
by  the  subscript  •. 

Since  localized  areas  of  interest  such  as  boundary  layers,  shear 
layers,  and  shocks,  are  present  in  the  flow  field,  a  non-uniform- 
grid  finite  difference  scheme  is  desirable.  This  will  allow  use  of  a 
fine  mesh  in  these  special  areas,  while  a  coarser  grid  can  be  used  in 
other  areas  preventing  machine  time  from  becoming  unwieldy. 
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A  second  order  accurate  centered  finite  difference  scheme 
for  a  variable  mesh  will  yield  approximations  for  both  the 
first  and  second  derivatives  which  depend  upon  the  values  of 
the  function  at  three  grid  points;  the  one  at  which  the  deriv¬ 
ative  is  required  and  the  two  adjacent  points.  If  such  a  scheme 
is  applied  to  the  steady  terms  of  the  continuity  equation,  for 
instance,  it  is  easily  shown  that  the  inclusion  of  the  central 
term  in  the  first  difference  approximation  yields  a  difference 
equation  which  is  not  in  conservation  form.  Therefore,  a  first 
order  accurate,  centered  difference  scheme,  which  yields  a  set 
of  difference  equations  in  conservation  form,  is  used.  For  two 
axial  steps  with  a  length  ratio  a  the  first  and  second  differences 
of  a  function  f  are 


n-1 

1 _ 

n 

J 

n+1 

1  - 

1 

Ax 

1 

a  Ax 

| 

**«>„ 


f  .  -  f  - 
n+1  n-1 

(a+1)  Ax 


<fxx>n 


fn*l  -  (1+a)fn+°  fn-l 
a (a+1) Ax8 


(2) 

(3) 


with  similar  expressions  for  the  radial  difference  with  the  radial 

step  size  ratio  /3  replacing  a. 

Taylor  series  expansions  of  f  ^  and  f^  about  f^  show  the 
error  of  the  above  schemes  with  respect  to  the  error  of  the  second 
order  central  scheme  for  a  uniform  mesh  to  be,  respectively 
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(4) 


V* 


and 


j  icrii  aa  n 

AX  (£  ) 

xxx  n 


4  (g-1)  ^xxx^n 

Lx  (£  ) 

xxxx  n 


Therefore,  the  accuracy  of  the  present  scheme  can  be  controlled 
by  keeping  the  relative  grid  spacing  a  close  to  unity. 

Since  variable  viscosity  and  heat  conductivity  are  being 
considered,  difference  approximations  for  terms  of  the  form 
(af  )  and  (af  )  will  be  required.  These  can  be  obtained  by 

XX  X  IT 

successive  applications  of  the  first  difference  expression  (2) 


^a^x^x^m,n  ^am,n+l  +  am,n^  *m,n+l 


-  [a  .  +  (1+a)  a  +  a a  ,  ]f 

m,n+i  m,n  m,n— l  m,n 


+  a  (a  +  a 
m,n  m 


(5) 


A  similar  expression  for  (af^)^  can  be  obtained  by  replacing  a  by 
B,  Ax  by  Ar,  and  Interchanging  m  and  n,  where  p  Is  the  ratio  of 
successive  mesh  Intervals  in  the  radial  direction.  The  cross 
derivative  expressions  are 

■  ■  m+1 
BAr 

-•  m 
Ar 
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^a^x^r^m,n  ^am+l,n^m+l,n+l  *m+l,n-l^ 


-  a  .  (f  .  ,.-  f  .  .)]/  (a+1)  (?+l)AxAr  (6) 

m-l,n  m-l,n+l  m-l,n-l  J/ 


r  (af  )  ]  =  [a  ^  (f  -  f  .  Al) 

r  x  m,n  L  m,n+l  m+l,n+l  m-l,n+l 


-  am,n-l(fm+l,n-l  "  f«-l.  n-l>  ]/ (o+1)  (P+1) *X*r 


(7) 


Time  Differencing 

The  form  of  the  spatial  difference  influences  the  choice  of 

a  time  differencing  scheme  through  stability  considerations.  One 

stable  scheme  consistent  with  the  above  spatial  differencing  is 

the  two-step  scheme  of  Brailovskaya.  There  are  two  stability 

criteria  on  time  step  associated  with  this  scheme,  one  linearly 

dependent  upon  the  spatial  step  size  (At  s  and  the  other  pro- 

c 

portional  to  the  square  of  the  spatial  step  size  divided  by  the 
kinematic  viscosity  (At  *  ) .  In  the  present  problem,  there  is 

a  large  expansion  of  the  flow  around  the  base  and  the  density  in 
the  base  region  becomes  small.  This  will  result  in  a  large  kine¬ 
matic  viscosity  and  severely  limits  the  time  step  results. 

4 

Allen  has  considered  a  modification  of  the  Brailovskaya 
scheme  which  removes  the  diffusion  dependent  stability  criteria. 

The  modification  consists  of  writing  the  independent  variables 
appearing  in  central  terms  of  twice  spatially  differenced  diffusive 
terms  at  the  updated  time  level  and  combining  them  with  the  time 
differenced  unsteady  terms  to  solve  for  the  independent  variables 
at  the  new  time  level.  This  works  because  we  are  only  interested  in 
the  steady  state  solution  and  how  we  get  there,  i.e.,  the  time 
evolution  of  the  solution  is  irrelevant. 
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Then  the  Wime  difference  scheme  for  the  present  differential 
equations  can  be  written  in  two  steps: 


1) 


p!+*  ■  pj!  _  +  Fi 


m,n  m,n  r  1 
in 


m,n 


m,n 


[(ou)^  +  ^  F  (p^u^v^e^u^X^)  ] 

m#n  r  z  m#n 

in 

jj  +  l  +  g  ,  j, 

m,n  r  Re  2  m,n 
m  ® 


C(pv)^l  +  ^  F  (D^u^v^,e^,u^,X^)  ] 

...  m,n  r  j  m,n 

vJ+i  ,  - SI - 

m,n 


v**l,  *  G  (u1) 

m,  n  it  K©  j  m,n 
m  ® 


(8) 


-j+1 

em#n 


{[p(e+q)]^  “<oq)^+  F4(pj,  uj,vj,e j+1) 
* _ m _ 


,i+1  ♦ 


^  G„  (kj) 


m,n  P  Re  r  4 
r  •  m 


m,n 


2) 


j+1  j  ,  A!  „  ,.j+l  -j+1  -j+1. 

^  n  *  n  r  P1  l*  »a  'V  L  n 

m#n  in  #  n  r  x  in  #  n 

in 


u 


j+i 

"m#n 


C (.«)!  •  *  *  »a(»,*l.«J*l.*Jrt.»,+1.uJrt.* j+1  >m  J 

m,n  r  z  m#n 

in 

C +  fsr 

m  ® 


,J+1 

m,n 


r/  J  ■  At  ^  /-j+1  „j+l  wj+l  -j+1  -j+1  z  j+l\  i 

r  (pv)^  _  +  tT"  F,  (PJ  ,eJ  ,UJ  ,X  J  )  ] 

hi  #  ri 


m#n  r  3 

_ BL 


0j+1  +  — G  (uj+1) 
m#n  r  Re  3  m,n 

m  • 


(9) 


m,i 
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j+1  l C»-  (e4<i)  ^  „} 


m,n 


j+i  +  -Is_  ii  G  (Ej+i, 

Dm,n  PrRe  r  4  'm.n 

r  os  m 


'here 


F  (pj,uj,vj) 
l  m,n 


1  * 

’  T77I777  r  <pru^m,n+l  '  <l>ru)m.n-l’ 
'  1771)77  C(orV)m+l.n  -  ^orv3m-l,n3 


P,(p3,u3,v^,e3,u3.X3) 

2  m,n 


“  TTT7 -  fpr(us+M  eXl]j  .  -  [pr  (u8  +M  eX)  ] j  J 

(a+l)Ax  VM  <»  Jm,n+1  *  Jm,n-lf 

-  Iff+ml  [(prUV)m+l,n  -  <»ruv)m-l,n3 


4/3  1  "j  1  1 

Re  1  a(a+l)Axa  ,n+l  m,n+l 


{ 


m^r 


(10) 


+  a[  (rn)^  n  1  n  1 

m,n  m,n-i  m,n-i 


[  +  (ru);J  „]  u; 


B(P+l)Ar2  l  L '~K,m+l,n  ~m+l,n 


+  B[(ru)^  +  (r|i)^  1  n^Um  1  n 

m,n  m- x  §  n  in™ x  $  n 
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2/3 


(a+1)  (3+1)  AxAr 


(v^  —  ) 

m,n+l  m+l,n+l  m-l,n+l 


-  (r*j)m,n-l(vm+l.n-l  '  Vm-l,n-l)  3 


-  "  -  - . -  f  (rii)  ^  ) 

(a+1) (P+1) AxAr  L  ^m+l#n  m+l,n+l  m+l,n-l 


-  (rn)-j  ,  (v^ 

m— l , n  m—  i , 


—  )  ] 
n+1  m-l,n-l; J 


"  (a+1) Ax  t(uv)m,n+l  “  (uV)m,n-l^  i 


F-i  (p^.u^v^e^n^X^)  _ 

j  m,n 


‘  (9+1) Ar"  ■f[‘>r(v>+“.eX)]m+l.n-  lortv’+iSeX)^)  } 


taTifc  C('iruv)m,n+1  -  (oruv)2,n-^ 


it  { rtsiiLr’  l - +  (ru)-  vJ 


m+l#n 


m,nJ  m+l,n 


+  fl[(ru)^  _  +  (ry)^  i  n3 

rn  /  n  nt*i  f  n  m~  x  f  n 


+  a(a+l)Axa  ^r^m,n+l  +  ^ru^m,n^  Vm, 


n+1 


+  a[ (ru)^  +  (ry)^  . ] 

m#n  nil  ri“i  nif  n*x 


(10  cont 
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2/3 


(a+1)  (8+l)AxAr  ^  ^ru^m+l,n  ^Um+l,n+l  Um+l,n-l^ 
-  (rw^m-l,n (um-l,n+l  “  Um-l.n-l}  3 


(a+1)  (fl+l)  Ax^r  ^r^m,n+l  ^Um+l,n+l  um-l,n+l^ 


-  (rM) .  ,  -  U-*  )] 

m,n-l  m+1,  n-1  m— i#n— l 


•  |)>mr  <“W>2-1.„> 


(2/3  )M 


m, 


-  (v* 


-  V* 


(3+1) ftr  m+l,n  m-l,n 


) 


(2/3  )M 


ULdL  (UJ 


-  u 


(a+1) ax  m,n+l  ra,n-l 

A  (mv)J  „  .  A  4 

4  m.n  /  ..i  1 

— - *—  >  +  (pM  eX) 

3  r  J  •  m,n 


) 


m 


p  (J  J  J  J  >3iu3+l  -j+1 

F4(p  ,u  ,v  ,e  ,u  ,k  / C*Y' u  'v  'm#n  “ 


fa? ITS  {Coru(\e+q)  ]^n+l  '  [cr^e+q>  } 


'  {[prV<>e+q)1m+l,n  '  Cpr^ye+q)]^!^  ] 


(10  cont'd) 
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+  oT  {«  (a+l>  **■  + 

i  00 


+  a[(rk)^  n  +  (rk);j  .]  (Ce)^  ] 

m,n  m,n-l  m»n-ii 


j 


B(6+l)&r*  |C(rk)m+l.n  +  <rk)2,n]  "e)m+l.n 


+  B[(rk>;>  „  +  (rk)3  3  (C®)^  ,  -  '!■ 

m,n  m-  J.  §  n  m*  x  $  n 


R ft  (t(tu)m,n+l+  (rw)m,n](u2,n+l)a 

00  * 

-  t(r“>m.n+l+(1+a» 


+■  rt[(ru)^  +  (ru)^  ,](u-J  ,)! 

m,n  m,n-±  m,n—i 


1/1 


j 


a(a+l)Bx>  |^>;,nn  +  (ru,ln]  lvln+l>3 


(10  cont’ 


-  C<r“>».n+l+(1+0> 


-  - 2/2 -  r(ruu)^  (vj  -  vj  ) 

(a+1) (fi+1) A*ArL v  ^  ;m,n+l  v  m+l,n+l  m-l,n+l; 

(ruvi)m^n_i  ^vm+i#n_l  Vm-l,n-l^ 


(a+1)  (p+1)axAx  ^rwV^m,n+l  ^Um+l,n+l  Um-l,n+l^ 


-  (ruv)m.n-l  <Um+l,n-l  ‘  "m-l.n-1^ 
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♦  [<uUV)m.n+l  ‘  (uUV>l.n-l] 

♦  7»wS  !c  (rw)^"  +  (ru,”-1-”]  <V"+1-n>S 

-  C(ru)i+l.n+(1+8)  (rM,i,n+  P(ru)m-l,n](Vn) 

+  8[<rw)^n  +  (r«>Ll.n](vm-l.n)S} 

♦  9  fc  (t(ru)m+l.n+  (r»)m.n^um+l.n>> 

-C<^>L,n+  d«)(ru)^n+  Mru,8.1>n]<i)S 

+  *tOw)i.„  +  (ru,»-l.»]  )  (10  concld) 


"  (a+l(  (0+D  Ax  A r  C  (ruv)m+l,n  (  m+l,n+l 


-  <*>*v,i-l.n  <um-l.n+l 


um+l,n-l> 


-  U  .  ,  J 


+  ■(^l)(^l)Mcirt(ruulm^.n  (Vm+l.n+l  '  viU.n-l’ 


-  ^r^u^m-l*n  ^vm-l*n+l  m-l#n-l 


Uuv’)*  n  -  (»*v)i-l.n3  > 


(B+l)  &r 


m+l#n 


(a+1)  AX*  C  (ru)m.n+l 


+  • (1+a)  (ru) 


^  +  a(ru) 


m,n-l' 


TR-763 


14 


Mife-  C<ru)>+li„  *  d«)(ru)^n  ♦  6  (ru)^.1>n] 


G,(uJ)_  n 

i  m,  n 


M0+lfc»  r<ru,m+l,n  +  (1+S)(r< r  +  e<r“>m-l.»’ 
a  (a+1)  AX>  C(«u)l,ntl  +  <^>«ru)^n  ♦  atr*)^] 


G  (kJ) 

4  m,n 


alaniTx’  [(rk)m,n+l  +  tl+tt)  <rk)m,n  +  o(rk)m.n-l] 
+  M«+ m?  C(rk,m+l.n  +  (l+B)  <rk,m.n  +  B(rk)m-l,n’ 


(11  concld) 


The  normalized  viscosities,  u,  are  obtained  from  e  through 

m,  n 


Sutherland's  law 


m,n 


(e  J 

m,n 


3/2 


(1+a)  /  (e  +a) 
m,n 


where  a  is  a  constant  and  the  normalized  thermal  conductivity,  k, 
from  the  definition  of  Pr  (for  a  perfect  gas  k  -  /t) . 

Note  that  the  same  functions  F  and  G  are  used  in  both  steps  of 

the  iteration  procedure  and  the  same  mesh  is  used  for  both  steps. 

1  4  5 

As  pointed  out  by  Cheng  and  Allen,  '  this  simplifies  the  program 
considerably. 
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JIM#**** 


In  the  actual  computer  programming  of  these  difference  equations 
the  F  functions  are  broken  into  the  sum  of  two  functions,  F  and  H. 

The  F  terms  contain  the  non-central  portions  of  the  twice  spatially 
differenced  diffusive  terms  whose  central  portions  have  been  used  at 
the  updated  time  level.  These  updated  central  terms  have  been 
transposed,  combined  with  the  left  hand  side  and  divided  through, 
appearing  as  the  G  functions.  The  F4  function,  which  is  part  of 
the  energy  equation,  is  somewhat  special  in  that  it  also  contains 

some  complete  differenced  velocity  terms  in  which  the  central 
terms  are  at  the  updated  time  level.  These  updated  velocities 
are  known  from  the  solution  of  the  momentum  equations. 

The  H  functions  contain  all  of  the  terms  which  are  not 
twice  differenced  diffusive  in  nature  and  therefore  have  no 
portions  at  the  updated  time  level.  These  terms  are  calculated 
by  general  differencing  subroutines,  which  use  difference  forms 
given  by  Eqs .  (2),  (6)  and  (  7  )  at  time  level  j. 

This  break  up  of  the  F  functions  removes  the  requirement 
of  writing  out  and  programming  the  difference  form  of  the 
complete  set  of  differential  equations.  Only  the  F  and  G 
functions  need  be  written  out  explicitly  while  the  H  terms 
are  calculated  from  the  general  subroutines. 

The  requirement  used  to  check  the  approach  of  the  un¬ 
steady  calculation  to  the  steady  state  is  that  the  density 
change  between  two  successive  time  steps  at  the  grid  point 
of  maximum  change  should  be  less  than  a  prescribed  small 
number . 


max 

m,n 


<  e 


(12) 
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BOUNDARIES 


The  preceding  set  of  difference  equations  is  used  to  compute 
the  flow  field  in  the  interior  of  the  domain  ABCFGHJA  of  Figure  1. 
It  remains  to  develop  procedures  for  treating  the  six  boundaries 
of  this  domain. 

Upstream 

As  noted  earler,  the  flow  field  along  the  upstream  boundary 
ABC  is  known  from  a  patching  of  an  inviscid  characteristics  solu¬ 
tion  and  a  boundary  layer  solution.^  The  values  of  the  independent 
variables;  p,  u#  v,  e,  from  these  solutions  are  used  as  fixed 
(with  time)  upstream  boundary  conditions  alcng  line  AC. 

Body  Wall  and  Base 

The  body  wall  and  base  boundaries,  AJ  and  JH  respectively,  can 
be  made  equivalent  to  each  other  by  interchanging  x  and  r.  In 
accordance  with  Allen's  numerical  results  of  negative  densities  on 
the  body,4  mesh  points  on  the  body  base  and  wall  are  being  avoided 
and  these  boundaries  are  taken  at  mid-mesh  locations  as  shown  below 
for  the  body  wall. 


BL  r 


4 r 
2 


Y  .  w 

7  /  /  /  rv  777  /  /  //  7  /  / 
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— ^Trn-"TTiTiiair«rr 


m**VK 


Alien's  emphasis  on  the  preservation  of  conservation  in  the  differ¬ 
ence  formulation  of  boundary  conditions  is  being  observed;  however, 
our  approach  to  this  differs  from  his  (which  appears  to  break  down 
for  a  nonuniform  mesh) .  We  accomplish  conservation  through  an 
inverse  consideration  of  telescoping  of  terms  under  quadrature 
approximation  to  integration.  That  is,  we  form  the  boundary  condi¬ 
tions  so  that  telescoping  of  internal  contributions  and  proper 
appearance  of  wall  terms  occurs  under  quadrature.  This  procedure 
should  automatically  insure  the  global  conservation  of  mass,  momentum 
and  energy. 

An  open  ended  trapezoidal  integration  of  the  difference  form 
is  carried  out  between  the  wall  and  the  outer  boundary,  i.e.,  for 
the  first  radial  difference: 


f  *  f  dr  -  f 

j  r  a 


f 

w 


_  f  j  f&r+ft£r  ,  fflAr+6Ar 

r'o1  2  1  r  1 1 '  2  1 


+  . . . 


(13) 


Point  1  and  those  above  it  are  general  points,  so  a  general  radial 
first  difference  form  similar  to  Eq.  (2)  is  used  for  f  L  f  I.  etc. 
Then  satisfaction  of  the  telescoping  requirement  near  the  wall  in 
Eq.  (13)  leads  to  a  three  point  first  difference  form  at  point  o 


f.  +f  -  2f 
_  I  _  l,n  o,n _ w,n 

r'o,n  (l+/3)Ar 


(14) 


Actually,  this  procedure  will  also  yield  a  special  difference  form 
for  f  I  ,  but  it  will  never  be  used  since  no  difference  equations 

f  09 

are  applied  at  the  cuter  boundary.  This  same  procedure  can  be  used 
for  the  second  derivatives  and  cross  derivatives. 
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(af)l  n  =  { CP-+2/9)  (a  +a  )  +  2a  ]  f 

r  r  *o,n  I  o , n  l , r  w , n  l.n 


(af  )  I 
x  r'o,n 


-  (l+2fl)  [a  +a  +  2  (1+2/5)  a  ]  f 

o  $  n  l  ,n  w  $  n  o  $  n 

+  8/90+1)  a  fw  n  }  /  /5(l+/5)  (1+2/3)  A ra 

( 1 5 ) 

^al,n^l,n+l  ^l,n-l^  +  ao,n^fo,n+l  fo,n-l^ 

-  2  aw,n(fw,n+rfw,n-l^/(1+a)(1+ft) 


(16) 


(af  )  I  =  [a  (f  Al+f.  ^  -  2  f  ) 

r  x'o.n  o,n+l  o,n+l  l,n+l  w,n+l 


“  *  n  l(f  n  l+fl  n  1  “  2  f „  „  n)]/ (1+00(1+0)  ****  (17) 

o,n-l  o#n-l  l,n-l  w,n-l 

A  possibly  more  accurate  scheme  can  be  obtained  by  using 
closed  ended  integration 


J*  f*dr  =  f  -  f  =  f  |  ^-  +  f  |  (^-  +  &£) 

J  <»  w  r'w  4  r'o  4  2 


w 


+  fr|i(fi£±4M)  +  ... 


(18) 


Point  1  and  those  above  it  are  still  general.  However,  a  second 
order  accurate  difference  formula  can  be  used  at  point  o 


-4 0*  f  +  (48s -1)  f  +f, 
f  |  =  _ _ _ S£d3 _ o>n  l,n 

r’o.n  B( 2/9+1)  Ar 


(19) 
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and  conservation  form  can  still  be  retained  by  solving  for  the 
f  lw  which  satisfies  Eq.  (18) . 

f  .  -W1-*)  fw.n+(4<’i-28-1)fo.J 

r'w,n  /JAr 


This  backed  out  difference  form  at  the  wall  is  never  actually 
used  since  no  difference  equations  are  used  along  the  wall  but 
a  half  mesh  out.  The  difference  forms  obtained  from  this  closed 
integration  scheme  were  programmed  but  no  calculations  were 
carried  out.  All  of  the  results  discussed  later  use  the  open 
ended  integration  difference  forms  given  by  Eq.  (14)  -  (17) . 

For  reference,  the  second  and  cross  difference  forms  for  the 
closed  integration  scheme  are  presented  here,  even  though  they 
have  not  yet  been  used  in  the  numerical  examples. 


(af  )  |  ={[(4tf-l)  a  +  a  ]  f 

r  r'o,n  l  o,n  l,n  l,n 


+  [  (1+2#?)  (4/1®  -6/3+1)  a  -a  -  S/^a  j  f 

o,n  1  § n  w f  n  o § n 

+  8fla[(l-«)  a  +fla  j  f  )  /  fP  (1+20)  Ara 
M  u  o,n  w,n  w,nj  K 

(af  )  I  =  la  A1[-4/ff  .  +  (4/3*  -1)  f  .,+f.  . .  ] 

r  x'o,n  I  o,n+lu  y  w,n+l  K  o,n+l  l,n+l 

-  a  .  C-4/58  f  n  1  +  (4/3a-I)  f  n  .  +  f.  n  ,]//5(l+a)  d+2/9)  AxAr 
o,n-l  w,n-l  o,n-l  l,n-l 
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(af  )  I  =  (-4^3  (f  -f  1)  +  (4/3a-l)  a  (f  -f  J 
x  r'o,n  l  w,n  w,n+l  w,n-l  H  o,n  o,n+l  o,n-l 


+  '  "d+at  (1+29)  ^  ■ 


Since  the  finite  difference  equations  are  not  applied  along  the 
wall  it  is  necessary  to  obtain  wall  values  by  other  means.  The 
velocity  components  u  and  v  are  both  known  to  be  zero  on  the  body 
wall  and  base.  For  other  quantities  such  as  p,  rp,  rk,  X  and  (, 
a  two  point  extrapolation  is  used 

f  =  [  ( 1+2/3)  f  -  f.  ]/20  (20) 

w,n  o,n  l,n 

The  wall  density  need  only  be  computed  at  the  final  time  if  the 

products  (pu)  and  (pv)  are  recognized  as  being  zero, 
w  t  n  w  f  n 

Either  of  two  wall  boundary  conditions  are  employed  to  obtain 

e  .  A  constant  wall  temperature  condition  can  be  applj  .  by  simply 
w  f  n 

setting 

e  =  constant  (21) 

w  i  n 

It  appears  that  it  would  be  possible  for  this  constant  to  vary  along 
the  wall,  but  this  was  not  attempted  here  for  the  sake  of  simplicity. 
For  an  adiabatic  wall  the  condition  that  de/dr  =  0  can  be  satisfied 
by  requiring 


e 

w,n 


e 


o,n 


(22) 


Both  of  these  options  are  provided  in  the  computer  program. 
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The  general  point  time  difference  scheme,  namely  Eqs.  (8)  and 

(9)  can  still  be  used  for  points  one  half  mesh  from  the  body.  The 

F  functions  are  broken  up  into  F  and  H  functions  as  described 

before.  The  H  functions  remain  the  same  as  the  general  case  except 

that  Eqs. (14)  to  (17)  are  used  to  compute  all  radial  differences. 

The  F  and  G  functions  for  points  one  haj.r  mesh  from  the  body  surface 

are  given  for  reference  in  Appendix  I.  The  quantities  (rji)  , 

w  f  n 

(rk)  and  P  are  obtained  from  the  extrapolation  formula  (20) 
w,n  w,n 

at  each  time  step.  The  difference  equations  along  a  line  of  gr^'d 
points  one  half  a  mesh  behind  the  body  base  can  be  similarly  ob¬ 
tained  or  can  be  evolved  from  the  body  wall  equations  by  inter¬ 
changing  x  and  r,  a  and  fi,  and  u  and  v. 

There  are  three  points  at  the  corner  formed  by  the  body  wall 
and  base,  which  require  special  considerations.  These  points  are 
designated  in  Figure  2  as  m,  n-1;  m,n  and  m-l,n.  They  are  special 
because  mixed  derivatives  require  difference  forms  of  both  side 
wall  and  base  wall  types.  Judicious  use  of  the  individual  differ¬ 
ence  forms  results  in  the  following  results  for  the  mixed  derivative 
differences  at  these  points. 

(af  )  I _ =  {a  (f  .  -f  .  ) 

r  x  m,nrl  L  m,n  m**!  ^ n 


-C(l+A)/fl]a  n  ,(f  ..  „  ,-f  n  ,)V (1+a)  (1+#) AxAr 
i  m,n-2  m+i,n-2  m,n-2  ) 


(23a) 
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^a^x^r^m,n-l  ^am+l#n-l  ^m+l,n  fm+l,n-2^ 


-  a_  _  ,  (f  _-f  _  2)]/U+a)/3  Ax^r 

m,n-i  m,n  nil  n™* 


^a^r^x^m,n  {am,n+l  ^m+l,n+l  ^m-l,n+l^ 


(23b) 


-[(!+*)/*]  Vn_1(fmtl(n.1-Vn-1)}/(^>M  AxAr 

«3c) 

^a^x^rlm,n~  {am+l,n  ^fm+l,n+l  ^m+l,n-l^ 

+[(l+a)/a]  am_1  n  (fm_l,n+l-fm-l,n)}/(1+a)(1+^)  *xAr 

(23d) 

^a^r^x^m-l#n  ~  ^am-l,n+l ^fm,n+l  fm-2,n+l^ 


-  a  If  -f  0  )  ]/a  (1+0)  Axa r 

m- 1 0  n  in  §  n  m—2  §  n 


(23e) 


(af  )  I  .  =  (a  (f  -f  .) 

x  r'm-l,n  l  m,n  m,n+l  m,n-l 


-  [  (l+a)  /a]  »n.2.n(fm-2,n+r£I»-2,n)}/(1+a)(1+',>iX4r 


(23  f ) 
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These  derivatives  only  appear  in  the  H  functions,  therefore  the 
appropriate  F  and  G  functions  can  be  used  at  each  point  with 
special  subroutines  for  the  mixed  derivatives  in  the  H  functions 
only. 


Centerline 

Along  the  boundary  represented  by  the  axis  of  symmetry  in  the 
wake, mesh  points  are  again  avoided.  A  special  finite  difference 
equation  for  use  one  mesh  away  from  the  axis  has  been  developed 
using  an  inverse  telescoping  procedure  similar  to  that  used  for 
the  wall  boundaries.  This  case  is  more  complex  since  the  radial 
difference  schemes  depend  on  whether  the  function  being  differenced 
is  odd  or  even  with  respect  to  radius  near  the  axis.  The  two  odd 
variables  are  r  and  v,  while  the  other  variables  p,u,e,p,etc. 
are  even.  Combinations  of  these  variables  are  considered  even  if 
they  contain  none  or  an  even  number  of  odd  variables  and  are  con¬ 
sidered  odd  if  they  contain  an  odd  number  of  odd  variables. 

By  definition,  an  even  function  has  zero  first  derivative  at 
the  axis,  r  =  0.  Therefore,  if  an  open  ended  integration  of  the 
first  difference  similar  to  Eq.  (13)  is  carried  out  between  the 
center  line  and  the  outer  boundary 
r 


f  f  dr  =  f  -f 
j  r  «b  c 


.(Ar+/g&r)  +  f  i  (fM± + 
1  2  r '  2 '  2  ' 


(24) 


and  th»  general  first  radial  difference  form,  equivalent  to  Eq.  (2) , 
is  used  for  points  2  and  above,  the  satisfaction  of  the  telescoping 
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requirement  in  Eq.  (24)  yields  a  special  first  difference  form 
for  point  1. 

f  L  „  "  <fi  +f,  -2f„  ) / (1+/3)  Ar  (f  even)  (25) 

it  If n  1  f n  2  f n  c § n 

For  an  odd  function  f  =  0  instead  of  f  |  .  Then  the  general  fcrm 
can  be  used  at  point  1  and  a  special  form  for  ^rlc  determined  by 
the  telescoping  requirement.  However,  since  no  difference  equations 
are  used  at  the  center  line  this  special  difference  form  is  never 
used.  Since  f  -  0,  the  first  difference  form  at  point  1  reduces  to 

fJi  „  -  fo  /<1+0)  *r  (f  odd)  (2S) 

r  i,n  4,n 

For  an  even  function  the  value  of  the  function  on  the  center 

line,  f  ,  is  required.  For  this  scheme  f  is  obtained  by  a  two 
c  c 

point  extrapolation 

f  Cfi  J1+0)B  -  f 9  2)  (f  even)  (27) 

Cf  n  x  f  n  &  f  n 
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The  higher  order  differences  can  be  obtained  in  a  similar 
manner.  Inspection  of  the  differential  equations  shows  that  all 
coefficients  a  in  derivatives  of  the  form  (a  fr)r  are  odd.  Using 
aca0,  the  telescoping  requirement  yields  the  second  difference 
forms 


(a£r’rll,n  *  'W  < V£l,/<,(8+1)  4r* 


(f  even) 


(28) 


(afr)r,l,n“{(ai+a2)f2"C(1+^)ai+a2;ifl}/fl(/5+1)Ar8 


(f  odd) 

(29) 


Eq.  (29)  is  the  general  form  with  f  *  0. 

The  cross  derivative  difference  forms  are  obtained  similarly, 
exercising  care  as  to  the  odd  or  even  nature  of  the  coefficients  a. 


^VxU.n  =  [al,n+l(£l,n+l+£2,n+l'2  ^.n+l’ 


-  ai,„-i(£i.„-i+£2.„-r2  W’/,ww  “4r 

(f  even)  (30) 


(a£r>xll.n-  Cal.n+l£2.n+ral.n-l£2.„-l’/(1+tt)(1+',)  4*4r 

(f  odd)  (31) 


(af  )  I  .  *  [a  (f  -f  ) 

x  r'l,n  l,n  l,n+l  l,n-l 


+  a2,n  (f2,n+l"f2,n-l) 


{ 


a  even,  f  even"t 
a  odd?r  f  odd  •» 


-  2a  if  ,-f  o  ) ]/ (1+a) (l+$)  Axa r 
c  i  n  c ,  n+i  c  f  n 


(32) 
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(aVrll,n  "  a2,n(f2,n+rf2,n-l)/(1+a<1+',)  ixar 

'a  even,  f  odd 


{ 


6r‘ - } 

a  odd,  f  evenJ 


(33) 


Eqs.  (25)  to  (33)  represent  the  difference  scheme  actually  pro¬ 
grammed  for  points  one  mesh  off  the  wake  center  line.  However, 
the  numerical  solution  showed  kinks  near  the  axes  in  the  radial 
density  plot.  The  extrapolation  for  center  line  values  of  even 
functions  was  judged  the  cause  of  these  kinks  and  a  new  center 
line  differencing  scheme  was  developed  but  is  not  incorporated 
into  the  computer  program  from  which  numerical  results  were 
obtained.  The  change  concerns  the  even  functions,  which  have 

fr  =0.  Instead  of  finding  f  by  telescoping  and  extrapolating 
c  r1 

for  fc#  the  general  difference  form  can  be  used  for  fr^  and  the 
telescoping  requirement  on  Eq.  (24)  used  to  find  fc.  This  results 
in 

f  =  f ^  (f  even)  (34) 


and  eliminates  the  need  for  extrapolations.  The  first  difference 
form  is  then 


£  I,  _  =  (f9  -f.  ) / (1+0)  Ar  (f  even) 

ri,n  z  ,n  i,n 


(35) 


The  second  difference  forms  are  unchanged  and  the  cross-derivatives 
which  are  changed  are 

(af  )  I.  =  [a  (f_  .,-f,  , , )  (f  even) 

r  x'l,n  u  l,n+l  2,n+l  l,n+l 

-  a.  n  .(f,  B  .-f.  b  jyil+aHl+o)  AxAr  (36) 
l,n-l  ^ ,n-i  l,n-i 


TR-763 


27 


(afx*rU.n  '  [a2,n(f2,n+l"f2,n-l) 


{ 


even,  f  odd 
6r 


a 

a  odd. 


f  even 


) 


-  a  (f  -f  . )  ]/(l+a)  (1+0)  AxAr  (37) 

c  §  n  X  §  n  *  i.  X  f  n  —  x 

At  the  final  time  step  Eq.  (27)  can  be  used  to  obtain  center  line 
values  of  the  even  functions.  The  F  and  G  functions  resulting 
from  using  these  difference  forms  in  the  differential  equations 
are  identical  for  either  scheme  and  only  the  subroutines  for 
computing  the  H  functions  are  affected. 

It  should  also  be  noted  that  a  completely  different  center 
line  differencing  scheme  can  be  obtained  by  use  of  closed  ended 
integration  in  the  telescoping  condition.  Most  of  these  improve¬ 
ments  were  programmed  in  preparation  for  a  trial  of  a  turbulent 
base  flow  calculation  (see  Appendix  m)  but  no  actual  calculations 
were  carried  out. 

The  F  and  G  functions  for  points  one  mesh  above  the  axis 
which  are  used  in  the  general  iteration  scheme,  are  given  in 
Appendix  I.  The  H  functions  are  the  same  as  a  general  point  with 
Eqs.  (25)  to  (33)  supplying  the  special  difference  subroutines  for 
the  center  line. 

Special  difference  formulas  are  needed  at  point  1,  0  (see 
Figure  3),  since  it  is  adjacent  to  both  the  base  wall  and  the  wake 
center  line.  The  first  and  second  differences  are  the  same  as  those 
along  the  boundaries  they  include.  The  axial  (x)  derivatives  are 
the  same  as  those  one  half  mesh  from  the  base  wall,  i.e.,  the  axial 
forms  of  Eqs.  (14)  and  (15),  while  the  radial  derivatives  are  the 
same  as  the  center  line  forms,  Eqs.  (25) -(27).  The  mixed  derivatives 
are  obtained  by  successive  applications  of  the  respective  first 
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differences.  These  are  listed  below  for  reference: 


(af  )  L  *  (a.  (f.  +f.  -2f  ) 

r  x'l,o  l  l,o  l,o  2,o  c,o 


+  iEll,l(fl,l+f2,l"2fc,l) 


(f  even) 

-  2#l.w'fl,w+f2.w-2fc.w>  iXW 


(38) 


“  {al,of2,o+al.lf2,l'2  al,wf2,w}/(1+o)(1+^  AxM 

(f  odd)  (39) 


(afx}rll,o  {al,o(fl,o+fl,r2fl,w) 

+  a2,o(f2,o+f2,l“2f2,w)  { 


a  even,  f  even., 
or 

a  odd,  f  odd 


} 


2ac,o^fc,o+fc,l“2fc,w)}/(l+a)  (1+0)  AxAr 


(40) 


(afx>rll,o  *  a2,o(f2,o+f2,r2f2,w)/(1+a)(1+a)  AxAr 

(a  odd,  f  even)  (41) 


Values  of  functions  at  point  c,  o  are  obtained  by  the  normal 

center  line  extrapolation  formula,  Eq.  (27) ,  while  functions  at 

point  1,  w  are  obtained  by  extrapolation  from  the  axial  form  of 

Eq.  (20).  All  the  functions  required  at  point  c,  w  contain 

u  or  v  which  are  zero  making  extrapolation  unnecessary. 
c,w  c,w 
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However,  it  would  be  of  interest  to  have  values  for  o  and 

vc,w 

e  for  the  final  steady  state  solution.  Except  for  the 
c ,  w 

constant  wall  temperature  case  where  e  =  e  ,  these  variables 

c,w  w 

have  been  left  undetermined  due  to  the  ambiguity  of  the  extrapola¬ 
tion  procedures. 

The  general  forms  of  the  H  functions  are  used  with  the 
appropriate  difference  forms.  This  special  point  requires  its 
own  set  of  F  and  G  functions  for  use  in  the  general  iteration 
scheme.  These  are  listed  in  Appendix  I. 

Downstream 

The  downstream  boundary  conditions  have  been  formulated  in  a 

three-point,  variable  mesh  extrapolation  along  lines  of  constant 
4 

radius.  Allen  found  this  procedure  to  be  satisfactory  if  the 
boundary  is  chosen  far  enough  downstream  for  the  flow  to  have 
returned  to  almost  completely  supersonic  speeds  (for  use  of  the 
downstream  data  line  as  input  to  a  far  wake  program  the  flow  must 
be  completely  supersonic) .  Another  restriction  is  that  the  angle 
between  the  streamline  and  the  constant  radius  line  at  each  down¬ 
stream  boundary  point  must  be  less  than  the  local  Mach  angle. 
Otherwise  the  upstream  points  on  the  constant  radius  line  will  not 
lie  in  the  zone  of  influence  of  the  boundary  point. 

Figure  4  shows  the  nomenclature  for  the  downstream  extrapola¬ 
tion.  The  second  order  accurate,  three  point,  variable  mesh 
extrapolation  formula  for  point  m,N  is 


fm.H  '  C4W4V4Vl>  (4ViVl+a’tN-2)  fm.N-l 


-  *XN 1+^-2  )<AVWN-l+4XN-2>  fm.N-2 

+  ‘VVl'Wl1  fm.N-3^4XN-l4XN-2<4XN-l+6XN-2) 

(42) 
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This  relation  can  be  used  to  obtain  the  values  of  each  variable 
at  each  time  step  for  both  parts  of  the  two  step  iteration. 

Upper  Boundary 

It  is  desired  to  obtain  values  of  the  independent  variables 
at  outer  most  radial  grid  (CF  in  Figure  1)  by  an  extrapolation 
procedure.  The  location  of  this  upper  boundary  line  is  picked  in 
such  a  way  that  the  flow  adjacent  to  it  is  inviscid  and  super¬ 
sonic.  The  downstream  boundary,  EG,  is  picked  where  the  center 
line  flow  is  expected  to  be  supersonic.  Then  the  first  family 
characteristic  line,  BF,  is  constructed  from  point  B,  which  is  in 
the  supersonic,  inviscid  region  of  the  upstream  data  line  above 
the  boundary  layer.  Point  F  is  determined  by  the  intersection  of 
this  line  and  the  downstream  boundary  line  EG.  Then  line  CF  will 
lie  completely  above  characteristic  line  BF  and  the  flow  adjacent 
to  it  will  be  inviscid  since  the  viscous  effects  generated  in  the 
boundary  layer  portion  of  the  inflow,  AB,  cannot  propagate  across 

the  characteristic  line. 

4 

Allen  used  an  extrapolation  procedure  which  involves  an 
inviscid,  steady  state  simple  wave  characteristic  solution.  The 
present  problem  requires  a  more  complex  procedure  due  to  the  axi- 
symmetry  and  the  nonuniform  flow  external  to  CF,  which  is  obtained 
from  an  independent  inviscid  flow  field  calculation.  The  current 
method  involves  a  steady  state,  two  step  characteristics  extrapola 
tion.  The  use  of  a  steady  state  solution  means  the  numerical 
results  will  be  incorrect  during  their  time  evolution,  but  when 
a  steady  state  is  reached,  it  will  be  correct. 

Figure  5  shows  typical  grid  points  near  the  upper  boundary. 
The  line  designated  by  M  represents  line  CF  in  Figure  1.  All 
points  up  to  line  M-l  are  computed  by  the  internal  difference  equa 
tions  and  it  is  desired  to  extrapolate  this  solution  to  a  typical 
point  M,n. 
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Before  describing  the  technique,  some  consideration  should 
be  given  as  to  just  how  much  data  at  M,n  should  be  prescribed 
from  the  known  mviscid  solution  and  how  much  extrapolated  from 
the  interior.  For  the  situation  of  Figure  5, data  are  transmitted 
out  of  the  calculation  region  at  B  along  two  lines;  the  streamline 
and  the  first  family  characteristic  line,  while  the  second  family 
line  at  B  transmits  data  into  the  region.  Each  of  the  character¬ 
istics  carries  one  piece  of  data,  via  the  characteristic  equation, 
while  the  streamline  carries  two  pieces  via  the  isentropic  and 
constant  stagnation  enthalpy  conditions.  Therefore,  for  this  case 
one  piece  of  data  is  prescribed  from  the  external  flow,  namely  the 
flow  deflection  9,  and  three  variables  are  extrapolated.  The 
numerical  examples  considered  later  coincide  with  this  geometry. 

If  a  problem  with  a  different  characteristics-streamline  arrange¬ 
ment  were  considered  the  computer  program  would  require  modification. 

The  characteristics  iteration  is  started  by  knowing  from  the 
interior  calculations,  all  the  variables  along  the  line  of  constant 
radius  represented. by  M-l  and  inputing  the  flow  deflection,  9,  at 
each  point  along  the  M  line.  Then  for  the  first  iteration  step  the 
first  family  characteristic  can  be  constructed  from  A  to  intercept 
the  M  line  at  B  (see  Reference  9) .  The  x  location  of  point  B  can 
be  found  from 


XB  - 

XA  +  (rM-rM-l,/tan(®AV 

(43) 

where 

u 

< 

a> 

tan"1(VA/UA) 

(43a) 

the  Mach  angle 

**A 

sin  ^ (1/M  ) 

A 

(43b) 

and  the  Mach  number 

Ma  =  [2qA/y(y-l)  e ^ 

(43c) 
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The  value  of  0  can  be  found  by  interpolating  between  points  M, 

B 

n  and  M,n+1. 

Next  the  streamline  from  B  back  toward  A  can  be  constructed 
using  an  average  of  slopes  at  A  and  B.  This  locates  point  C. 

rc  ■  rM  -  (xB‘*n-l)  (tan  Vtan  V/2  <44> 

It  is  assumed  that  the  variables  at  point  M,  n-1  are  known  from 
the  preceding  extrapolation.  If  n-1  is  the  upstream  boundary 
they  would  be  prescribed.  Then  the  values  of  p,e  and  q  at  C  can 
be  found  by  interpolation  between  points  A  and  M,n»l. 

Two  conditions  are  known  to  apply  along  the  streamline  CB.  ^ 
The  first,  which  requires  the  flew  to  be  isentropic  along  the 
streamline  can  be  expressed  by  relating  the  normalized  pressure 
and  density  as  follows 

P  =  d  py  (45a) 

where  d  is  a  constant  and  the  equation  of  state  gives  P  =  X  pe. 
Elimination  of  P  allows  Eq.  (45a)  to  be  written 

d  =  X  e  (p) 1_y  (45b) 

The  second  condition  is  that  the  stagnation  enthalpy  is  constant 
along  a  streamline 

H  *  e  +  q  ■  constant  (46) 

The  values  of  the  constants  d  and  H  can  be  evaluated  by  applying 
equations  (45b)  and  (46)  at  point  C.  These  will  be  used  later  at 
point  B. 


TR-763 


33 


9 

Next  the  first  family  characteristic  equation  along  line 
AB  can  be  used  to  obtain  the  normalized  pressure  at  B.  For  the 
first  iteration  step  this  equation  is 

Pfi  ■  -  [2yaPA/(y+l)  ]  (  (i6b“6a)/cos  pA  sin 

+  [sin  8A/rA  008  MA  cos  (®A+MA)  3  (XB-XA) }  (47) 


Since  d  and  H  are  constant  along  streamline  CB,  their  values  at 

B  are  the  same  as  at  C  so  Eqs.  (45a),  (46)  and  the  state  equation 

can  be  used  at  B  to  determine  p  ,  e  and  a  from  P  . 

B  B  d  B 

pb  -  <Vd)1/v  (48) 


eB  *  PB/',BXB  (49) 

%  ~  H  -  ®B  (5°> 


These  quantities  together  with  80  represent  the  first  step  of  the 
iteration  solution  at  point  B.  For  the  second  step  the  process  is 
repeated  using  average  values  of  the  slopes  at  both  ends  of  the 
characteristic  and  stream  lines.  This  more  accurately  locates 
points  B  and  C.  Equation  (47)  is  recomputed  using  similarly 
averaged  coefficients  of  the  8  and  x  terms  and  final  values  of 

Pb' 

is  completed,  the  velocities  u  and  v  can  be  obtained  from 


e, 


B'  S 


and 


are  determined.  After  the  second  iteration  step 
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UB  =  (2  V"-’  H  cos  8b 

.  i.  <51> 

VB  =  U%/MJ  Sin  ®B 

The  values  of  the  four  independent  variables  p,u#v  and  e  can  be 
obtained  at  M,n  by  interpolating  between  points  M,n-1  and  B. 

In  an  actual  problem  point  B  may  not  fall  between  line  n 
and  n+1  but  may  be  between  n+1  and  n+2  or  between  n-1  and  n,  etc. 
The  computer  program  has  taken  these  possibilities  into  account 
by  using  various  interpolations  and  or  extrapolations.  The  only 
geometric  restriction  on  the  flow  is  that  one  characteristic 
line  and  the  streamline  lie  below  the  upper  boundary  line,  while 
the  other  characteristic  line  lies  above.  The  program  can  be 
adapted  to  other  cases  but  this  has  not  been  done  as  yet. 
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NUMERICAL  RESULTS 


Computer  Program 

The  finite  difference  scheme  described  in  the  preceding  sections 
was  programmed  in  FORTRAN  Version  2.3  (includes  all  features  of 
FORTRAN  IV  plus  some  additions)  for  the  CDC  6400/6500/6600  computers. 
Three  grid  configurations  were  used: 


Configuration 

1 

2 

3 

r  direction  grid  points 

10 

50 

50 

x  direction  grid  points 

20 

50 

70 

Total  points  in  mesh 

170 

2040 

3040 

Execution  time  for  2  half-steps 
on  CDC  6600 

1.4 

sec  15 

sec  20. ( 

The  total  number  of  points  for  any  configuration  is  equal  to  the 
product  of  the  numbers  of  points  in  the  x  and  r  directions  minus 
those  points  which  fall  within  the  body. 

It  is  our  estimate  that  the  maximum  size  grid  configuration 
that  could  be  run  (by  expanding  the  appropriate  program  dimensions) 
is  90  x  90  or  some  other  combination  of  x  and  r  grid  points  whose 
product  is  equivalent.  This  maximum  grid  size  estimate  is  based 
upon  the  program  volume,  plus  the  storage  necessary  for  the  com¬ 
puter  operating  system,  plus  the  additional  storage  needed  at  load 
time.  The  cost  of  running  the  expanded  grid  is  approximately  pre¬ 
dictable  from  the  above  table  by  equating  the  ratios  of  execution 
time  and  total  mesh  points.  There  is  another  cost  factor.  The 
number  of  iterations  to  reach  convergence  is  greater  with  the 
finer  mesh  because  a  disturbance  can  propagate  at  most  one  grid 
point  in  one  half  step. 
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A  separate  computer  program  was  written  to  take  data  from  the 
intermediate  data  files  (TAPE  7)  and  create  digital  plots  of 
velocity  vectors,  streamlines  and  pressure  contours  in  the  flow 
field.  The  plotting  war  performed  on  a  Calcomp  Plotter.  The 
CDC  6600  software  required  to  run  this  program  is  the  Calcomp 
Plotter  Software  package;  specifically,  subroutines  PL0T,  SYMBOL 
and  NUMBER.  Instructions  for  data  input  to  the  plotting  program 
are  given  in  Appendix  II... 

Input  Data 

Essentially,  the  required  input  consists  of  the  grid  geometry, 
flow  parameters  and  initial  data  for  the  entire  grid.  A  complete 
card  file  input  description  is  given  in  Appendix  II.. 

The  first  step  in  choosing  the  grid  configuration  is  the 
location  of  the  boundaries  of  the  numerical  calculations.  As  out¬ 
lined  in  the  section  on  boundaries,  the  upstream  boundary  is  placed 
about  five  boundary  layer  thicknesses  upstream  of  the  base,  the 
downstream  boundary  is  placed  at  a  location  where  the  flow  is 
ecpected  to  have  returned  to  completely  supersonic  and  the  upper 
boundary  is  located  so  that  it  is  in  a  completely  inviscid  flow 
regime.  In  picking  a  mesh  size  it  must  be  remembered  that  the 
object  is  to  compute  flow  fields  of  interest  with  as  much  accuracy 
as  possible.  Since  the  spatial  stability  is  dependent  upon  the 

Reynolds  number  per  mesh  (Re  ) ,  too  large  a  value  causes  the  program 

A 

to  fail.  Therefore,  the  use  of  a  coarse  mesh  requires  a  restriction 

of  the  magnitude  of  the  free  stream  Reynolds  number  (RS^) ,  in  order 

to  limit  Re  .  On  the  other  hand,  the  interesting  features  of  the 
A 

flow,  such  as  the  recirculation  region,  tend  to  be  stretched  out 
over  a  greater  area  as  Re  is  increased.  This  shows  the  opposing 

GD 

effects  of  trying  to  have  the  recirculation  region  large  enough  to 
be  observable,  which  requires  large  Re  ,  and  trying  to  use  a  coarse 

CD 
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mesh  to  cut  running  costs ,  which  requires  a  lower  Re^  to  give 
spatial  stability.  Furthermore,  the  grid  must  be  fine  enough 
to  resolve  the  resulting  flow  field  (e.g.,  the  gradients  in  the 
body  boundary  layer) .  The  axial  step  size  is  usually  taken 
larger  than  the  radial  since  smaller  gradients  are  expected  in 
the  axial  direction.  Care  must  be  exercised  in  using  the 
variable  mesh  feature  since  the  magnitude  of  s  or  ^  affects  the 
order  of  accuracy  of  the  difference  formulas  (Eq.  (4)).  The 
final  configurations  must  represent  a  compromise  between  these 
considerations . 

The  flow  conditions  are  described  by  specifying  the  free 
stream  (ahead  of  shock)  Mach  number,  Reynolds  number,  M  =ua/e  , 

00  00  CD 

V Pr<#  and  for  the  constant  wall  temperature  case,  e^/e*  (see 
Appendix  II) .  It  is  also  necessary  to  give  the  values  of  the 
four  independent  variables  on  the  upstream  boundary  and  the  flow 
deflection  along  the  upper  boundary. 

To  start  a  calculation,  initial  values  of  each  independent 

variable  must  be  known  for  each  grid  point  at  which  the  finite 

difference  equations  are  applied.  The  very  first  time  a  case  is 

run  a  guess  is  made  for  this  data.  This  can  be  based  on  the  given 

upstream  data,  a  known  inviscid  flow  field  or  can  be  a  uniform 

flow,  etc.  In  practice,  the  computation  is  run  for  a  set  number 

of  time  steps,  and  the  results  observed.  Then  another  series  of 

time  steps  is  run  using  the  results  of  the  previous  computation 

as  initial  data.  Means  of  doing  this  are  described  in  Appendix  II.. .. 

When  it  is  desired  to  vary  flow  parameters,  say  Re  ,  itis  advantageous 

00 

to  start  from  some  intermediate  step  in  a  previously  run  case  whose 
Reynolds  number  is  closest  to  the  desired.  Also,  if  a  finer  mesh 
is  desired,  the  output  from  a  coarse  mesh  run  can  be  interpolated 
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and  used  as  initial  data.  The  machine  output  at  the  steady 
state  condition  has  the  same  form  as  the  initial  data  as 
described  in  Appendix  II. 

Time  Step  and  Steady  State  Criteria 

The  time  step  size  is  related  to  the  spatial  mesh  size  by 
the  stability  condition.  Normally  there  are  two  stability  con¬ 
ditions,  howevei ,  it  is  assumed  that  the  differencing  scheme 
used  here  has  eliminated  the  viscosity  dependent  condition.  An 
approximation  to  the  other  stability  condition  can  be  obtained  by 

considering  the  differential  equations  with  inviscid  terms  only 

(4 ) 

and  linearizing  them.  Following  the  procedure  of  Allen  the 
axisymmetric  stability  condition  is 


At  *  Ar  /  {(££)  ^  |u|  +  |v| 


Equation  (52)  is  computed  for  each  point  at  each  step  of  the  itera¬ 
tion.  To  be  sure  of  stability,  the  program  uses  a  At  equal  to 
85%  of  the  minimum  value  obtained  from  Eq.  (52) .  For  the  coarse 
mesh  (170  points).  At  «  0.13  to  0.14,  while  for  the  fine  mesh 
cases  (2040  -  3040  points).  At  &  0.03  to  0.04. 

The  solution  is  assumed  to  have  reached  the  desired  steady 
state  when  the  density  changes  by  less  than  a  small  number  « 

between  successive  time  steps  as  given  by  Eq.  (12).  For  the  casejs 

-5 

considered  here,  «  **■  10  was  found  satisfactory. 
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Numerical  Examples 

The  computer  program  for  a  spherically  capped  cylindrical 
body  was  run  initially  with  a  coarse  mesh  grid  for  various 
Reynolds  number  flows.  Even  though  they  lack  accuracy  these 
coarse  mesh  runs,  which  are  relatively  inexpensive,  can  be  used 
to  learn  much  about  the  program's  limitations,  such  as  maximum 
Reynolds  number  for  a  given  mesh  size,  before  running  the  more 
expensive  fine  mesh  cases. 

All  the  cases  considered  had  a  free  stream  Mach  number, 

M  =12.5,  ufl/e  =78.8,  y  =  y  =1.3,  and  e  /e  =*6  for  constant 

•  00  00  00  W  ® 

wall  temperature  cases.  The  Reynolds  number,  Re^,  was  varied 
between  40  and  10,000  and  two  values  of  the  Prandtl  number  were 
used,  Pr  =  Pr  =  1.0  and  0.72.  An  independent  inviscid  character- 

OB 

istics  calculation  was  used  to  obtain  the  upper  boundary  flow 
deflection  and  the  inviscid  portion  of  the  upstream  boundary  line. 

A  boundary  layer  calculation  was  made  to  complete  the  upstream 
data  requirements. 

It  should  be  noted  that  all  dimensions  are  in  terms  of  the 
base  half  height  or  cylinder  radius,  and  *  refers  to  the  uniform 
conditions  outside  the  bow  shock.  The  lower  Reynolds  number  cases 
were  treated  with  a  coarse  grid  having  170  points,  while  a  fine 
mesh  having  2040  points  was  used  for  the  larger  Reynolds  numbers. 

Coarse  Mesh  Solutions  :  The  first  runs  were  made  with  a 
coarse  grid  for  the  constant  wall  temperature  conditions.  A  low  Rey¬ 
nolds  No. case.  Re  =  40  was  used  to  see  if  convergence  to  a  steady 

CD 

state  solution  was  obtainable.  A  steady  state  was  approached  for 
this  case  after  many  parameter  adjustments  required  in  getting 
familiar  with  the  program  operation.  However,  for  this  coarse 
mesh  the  recirculation  region  was  not  observable.  Therefore,  a 
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case  for  Rew  =  400  was  run.  A  converged  solution  was  obtained  but 

the  recirculation  region  was  still  not  evident.  A  run  was  then 

made  with  Re  =  1200.  This  solution  converged  and  showed  a 
00 

definite  recirculation  region.  This  is  evidence  of  the  mesh 
size  tradeoff  mentioned  previously.  If  the  Reynolds  number  is 
kept  small  to  insure  stability  for  a  coarse  mesh,  which  is  cheaper 
to  run,  the  recirculation  region  will  be  small  and  the  coarse  mesh 
will  not  intercept  it.  It  should  be  mentioned  that  some  of  these 
coarse  mesh  runs  have  small  negative  density  values  behind  the 
base.  These  seem  to  be  localized  and  do  not  affect  the  rest  of  the 
flow  field.  They  tend  to  disappear  as  steady  state  is  approached, 
especially  for  the  finer  mesh  cases. 

Once  it  was  established  that  a  recirculation  region  could  be 
obtained,  the  next  step  was  to  increase  the  Reynolds  number  to 
find  the  limit  for  which  a  solution  could  be  obtained  for  the  coarse 
mesh.  After  an  Re  9  1800  case  was  run  to  convergence  a  value  of 

GD 

ReB  =  4000  was  attempted.  For  this  later  case  it  was  not  possible 
to  obtain  a  steady  state  condition.  The  calculation  terminated 
itself  when  the  internal  energy  became  negative  in  a  small  region 
of  the  grid.  This  is  a  positive  definite  quantity  for  which  negative 
values  are  unacceptable.  All  computations  for  which  steady  state 
solutions  could  not  be  reached  were  terminated  in  this  manner.  Sub¬ 
sequently,  a  case  with  Re  =  2500  was  run.  This  case  also  failed, 

QD 

giving  a  limiting  Re  of  between  1800  and  2500  for  coarse  grid. 

Since  the  terminating  factor  appears  to  be  negative  internal 
energy  it  would  seem  likely  that  the  upper  limit  of  Rew  could  be 
increased  by  increasing  the  effect  of  the  energy  dissipation  terms. 
This  was  investigated  by  reducing  the  Prandtl  number,  which  appears 

I 

in  the  denominator  of  these  terms,  from  1.0  to  0.72.  When  the 
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Re^  =  2500  case  was  rerun  for  Pr  =  0.72,  a  converged  solution 

was  obtained.  This  was  also  true  for  Re  *  3000;  however,  an 

00 

Re^  =  3500  case  failed.  Therefore,  the  Reynolds  number  limit  for 
a  given  mesh  can  be  increased  by  enhancing  the  effect  of  the 
energy  dissipation  terms.  All  subsequent  calculations  were  run 
with  Pr  »  0.72.  Many  additional  partial  runs  with  different 
Reynolds  numbers  and  coarse  grid  arrangements  were  initiated  in 
the  course  of  the  learning  process.  These  are  not  reported  here 
since  they  do  not  lead  directly  to  any  conclusions. 

Fine  Mesh  Solutions:  The  coarse  grid  cases  discussed  above 
cannot  be  expected  to  give  accurate  results.  They  have  been  used 
to  investigate  the  feasibility  of  obtaining  solutions  and  the 
nature  of  the  breakdown  at  large  Reynolds  numbers.  In  order  to 
obtain  a  usable  solution  a  fine  mesh  of  50  x  50  points  was  con¬ 
sidered.  Due  to  the  body  location  there  are  actually  only  2040 
computation  points.  Running  times  to  reach  steady  state  for  this 
mesh  depend  upon  initial  data  but  were  between  1  and  1-1/2  hours 
for  the  cases  considered.  The  initial  case  run  on  this  mesh  was 
the  largest  Reynolds  number  case  obtained  for  the  coarse  mesh; 

Re^  ■  3000.  A  convergent  solution  was  obtained. 

An  analysis  of  Burgers'  equation  in  one-dimension  shows  that 
for  local  Reynolds  numbers,  based  on  step  size,  greater  than  2.0, 
the  numerical  solution  exhibts  point-to-point  oscillations. 
However,  for  our  more  complex  system  and  geometry,  this  does  not 
seem  to  be  true.  Point-to-point  oscillations  do  arise,  however, 
for  this  Re  =  3000  fine  mesh  case,  and  accompany  large  curvatures 
in  the  solution. 
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Figure  6  presents  radial  plots  of  the  density  at  stations 
upstream  of  the  corner.  The  solid  curve  without  circles  is  the 
upstream  starting  profile,  and  the  circles  represent  computed 
points  somewhat  downstream.  In  the  left  hand  curve  for  constant 
wall  temperature,  the  computed  points  show  point-to-point  oscilla¬ 
tions  near  the  body  surface  even  though  Re  is  smaller  there  than 

Ar 

at  larger  values  of  r.  However,  in  the  range  of  radius  r  where 

oscillations  occur,  the  starting  profile  and  the  computed  profiles 

have  large  curvatures.  It  was  therefore  concluded  that  point-to- 

point  oscillations  can  occur  when  the  mesh  is  not  fine  enough  to 

represent  curvature  or  local  changes  in  a  variable. 

To  check  this  conclusion,  an  adiabatic-wall  case,  which  has 

a  smoother  density  profile,  was  run.  The  right  hand  side  of  Figure  6 

shows  that  no  oscillations  are  present.  The  values  of  Re  in  this 

Ar 

case  were  small,  and  in  order  to  further  verify  the  dependence  of 
oscillations  on  large  curvatures  the  same  case  was  run  at  a  larger 
free-stream  Reynolds  number  (10,000  instead  of  3,000).  This  computation 
and  the  results  appear  in  Figure  7.  Now  Re^r  is  large  but  no  oscilla¬ 
tions  exist  since  the  curvatures  are  small. 

Figure  8  shows  the  axial  variation  of  density  just  above  the 
body  surface.  The  upper  curves  are  for  the  constant  wall  tempera¬ 
ture  case  and  show  oscillations,  evidently  remnants  of  the  radial 
oscillations.  The  two  curves  show  the  time  evolution  of  the 
solution.  The  lower  curves  are  for  the  adiabatic  wall,  Re^-lOjOOO 
case  and  show  only  slight  oscillations  even  though  values  of  Re^r 

are  about  the  same  as  for  the  constant  T  case. 

w 

Figure  9  shows  radial  plots  of  density  and  internal  energy 

for  the  constant  T  case  at  a  location  one  half  height  behind  the 

w  3 

base.  The  density  oscillations  present  in  the  body  boundary  layer 
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have  propagated  downstream.  The  adiabatic  wall.  Re  =  10,000 

99 

case  is  plotted  in  Figure  10  at  the  same  station  and  contains 
no  oscillations.  The  time  evolution  of  p  is  also  indicated. 

Figure  11  is  a  similar  plot  for  the  adiabatic  wall.  Re  =10,000 

CD 

at  a  station  one  diameter  behind  the  base. 

Figures  12,  13  and  14  are  additional  machine  plotted 
results  (see  Appendix  II)  of  the  adiabatic  wall.  Re  =10,000  case. 
Figure  12  gives  streamline  contours.  The  recirculation  region  is 
well  defined  by  the  closed  contours,  which  are  plotted  at  finer 
increments  of  the  stream  function  than  the  outer  flow.  The  stag¬ 
nation  point  is  indicated  by  the  intersection  of  the  zero  stream¬ 
line  and  the  flow  center  line  (r  =  0) . 

The  lines  in  Figure  13  represent  the  velocity  vectors  through¬ 
out  the  computed  field.  Each  vector  starts  at  the  small  cross, 
each  of  which  represents  a  grid  point.  The  stagnation  point  is 
represented  by  the  asterisk. 

Machine-plotted  pressure  contours  are  shown  in  Figure  14.  The 
convergence  of  the  contours  at  the  downstream  end  of  the  calcula¬ 
tion  field  evidently  represents  the  formation  of  a  recompression 
shock. 

The  increasing  of  Re  from  3,000  to  10,000  caused  the  Mach 
number  at  the  downstream  end  of  the  grid  to  become  subsonic. 

Hence  a  computation  with  a  longer  domain  (50  x  70)  was  carried 
out  to  insure  supersonic  flow  at  the  downstream  boundary.  The 
results  in  Figures  15  and  16  show  that  the  curves  plotted  for  the 
short-domain  calculation  are  only  slightly  affected. 
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SUMMARY 


GASL  has  developed  a  computer  program  which  treats  the 
near  wake  of  an  axisymmetric  body  having  a  truncated  cylindri¬ 
cal  base.  The  numerical  approach  is  a  two-step  time-dependent 
finite  difference  procedure  for  the  compressible  Navier  Stokes 
equations . 

The  pregram  has  the  following  features: 

1)  Axisymmetric  flow  field. 

2)  The  calculated  flow  field  includes  the 
region  upstream  of  the  base  where  corner 
effects  are  present. 

3)  The  inviscid  flow  along  the  outer  hori¬ 
zontal  boundary  may  be  nonuniform. 

4)  A  variable  mesh  is  used  to  allow  higher 
resolution  in  areas  of  flow  nonuniformities. 

5)  Variable  viscosity  and  conductivity  are 
admitted. 

6)  Options  of  adiabatic  wall  or  fixed  wall 
temperature  are  admitted. 

7)  Both  the  equation  system  and  the  finite 
differencing  are  done  in  conservation  form. 

8)  The  stability  condition  for  the  time-dependent 
solution  is  independent  of  the  viscosity 
(Cheng-Alien  scheme) . 

9)  Difference  forms  at  the  boundaries  are 
obtained  by  means  of  telescoping  with  the 
internal  difference  forms. 
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A  number  of  numerical  examples  have  been  run.  Calculations 
using  a  coarse  grid  have  shown  that  for  small  Reynolds  numbers, 
where  the  recirculation  region  is  small,  the  recirculation  region 
may  not  be  intercepted  by  the  grid.  For  larger  Reynolds  numbers 
the  recirculation  region  becomes  evident,  however,  a  limit  of 
Reynolds  number  is  reached  for  which  a  given  grid  size  will  not 
yield  a  convergent,  steady  state  solution.  This  upper  Reynolds 
number  limit  can  be  increased  by  enhancing  the  effect  of  the 
energy  dissipation;  by  decreasing  the  Prandtl  number,  for 
instance. 

Point-to-point  oscillations  in  fine  mesh  solutions  do  not 
necessarily  arise  when  the  local  Reynolds  number,  based  on  step 
size,  is  greater  than  2.0  as  is  the  cane  for  the  one-dimensional 
Burger's  equation.  However,  such  oscillations  do  arise  when  the 
mesh  is  not  fine  enough  to  accurately  represent  curvatures  in 
tho  solution  variables.  This  was  verified  by  obtaining  smooth 
solutions  for  an  adiabatic  wall  case,  with  small  variable  curva¬ 
tures,  for  either  the  same  free  stream  Reynolds  number  or  local 
Reynolds  number  per  mesh  as  a  constant  wall  temperature  case  with 
large  curvatures,  viiich  exhibited  oscillations. 
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